function test_example_testing_bem_created_leadfields

% DATA no
% MEM 2gb
% WALLTIME 00:10:00

%
%% Testing BEM created EEG lead fields
%
%% # Introduction
%
% The following section compares lead fields calculated by means of BEM implementation and theoretical model, by correlating their values for different densities of the head model vertices. It is expected that solutions for a single shell EEG realistic head model converge to the theoretical one with an increasing number of points in the mesh.
%
% For the simplest case, the BEM and the theoretical solutions for EEG lead field can be compared directly, because theoretical EEG forward solution for a spherical shell is known. Let's calculate the electric field on top of a sphere, generated by one dipole source with different positions along the z-axis:
%
% Create a spherical volume conductor of radius 100, conductivity 1 and center [0,0,0]
vol   = [];
vol.r = 100;
vol.c = 1;
vol.o = [0 0 0];
vol.unit = 'mm';

% Create a set of electrodes on the upper half of a sphere
[X, Y, Z] = sphere(10);
pos = unique([X(:) Y(:) Z(:)], 'rows');
pos = pos(pos(:,3)>=0,:);

elec = [];
elec.elecpos = vol.r * pos;
elec.label = {};
nelec = size(pos,1);
for ii = 1:nelec
  elec.label{ii} = sprintf('vertex%03d', ii);
end

% Define positions of dipoles (along z axis)
zp = linspace(0,100,50)';
pos = [zeros(size(zp)) zeros(size(zp)) zp];

% Define the corresponding spatial grid
cfg = [];
cfg.inwardshift = 5;
cfg.sourcemodel.pos = pos;
cfg.headmodel = vol;
cfg.elec = elec;
sourcemodel = ft_prepare_sourcemodel(cfg);

%
%% # Building the geometrical head model with BEM
%
% The following code generates different BEM (Dipoli) system matrixes, describing head properties which are depending on the increasing number of triangles used for each mesh. This implements Oostendorp and van Oosterom, 1989, eq. 2.3, 2.4
%
% construct a dense triangulated mesh
[X, Y, Z] = sphere(100);
pos = unique([X(:) Y(:) Z(:)], 'rows');
mesh.pos = pos*100;
mesh.tri = convhulln(mesh.pos);

vertic = [50 80 100 200 2000];
for ll=1:length(vertic)
  cfg = [];
  cfg.method = 'dipoli';
  cfg.conductivity  = 1;
  cfg.numvertices = vertic(ll);
  cfg.isolatedsource = false;
  volbem{ll} = ft_prepare_headmodel(cfg, mesh);
end

%% # Simulation
%
% The following steps generate the leadfields with
%
% # the BEM model integrated in the forward equations,
% # the theoretical forward solution for a dipole in spherical conductor.
%
% The dipole positions are defined in the variable |sourcemodel|, the headmodels are defined in `volbem{ll}` and the sensor positions in |elec|.
%
% calculate BEM leadfield
for ll=1:length(vertic)
  cfg = [];
  cfg.sourcemodel = sourcemodel;
  cfg.headmodel = volbem{ll};
  cfg.elec = elec;
  leadfieldBEM{ll} = ft_prepare_leadfield(cfg);
end

% calculate theoretical single sphere leadfield
cfg = [];
cfg.grid = sourcemodel;
cfg.headmodel = vol;
cfg.elec = elec;
leadfieldSphere = ft_prepare_leadfield(cfg);

%% # Result
%
%
% The pinky arrow describes the correlation curves of meshes with increasing number of triangles. The last mesh (2000 vertices) has a flat correlation curve at value y=1.
